#include "hip/hip_runtime.h"
#include <hip/hip_runtime.h>
#include <stdio.h>
#include "kernel.h"
#include "util.h"
#define TX 32
#define TY 32
#define RAD 1
#define TPB 64
#define ATOMIC 1
#define REAL double

//__global__
//void dotKernel(REAL *d_res, const REAL *d_a, const REAL *d_b, int n) {
//	const int idx = threadIdx.x + blockDim.x * blockIdx.x;
//	if (idx >= n) return;
//	const int s_idx = threadIdx.x;
//
//	__shared__ REAL s_prod[TPB];
//	s_prod[s_idx] = d_a[idx] * d_b[idx];
//	__syncthreads();
//
//	if (s_idx == 0) {
//		int blockSum = 0;
//		for (int j = 0; j < blockDim.x; ++j) {
//			blockSum += s_prod[j];
//		}
//		printf("1Block_%d, blockSum = %d\n", blockIdx.x, blockSum);
//		// Try each of two versions of adding to the accumulator
//		if (ATOMIC) {
//			atomicAdd(d_res, blockSum);
//		}
//		else {
//			*d_res += blockSum;
//		}
//	}
//}

//void dotLauncher(REAL *res, const REAL *a, const REAL *b, int n) 
//{
//	REAL *d_res;
//	REAL *d_a = 0;
//	REAL *d_b = 0;
//
//	hipMalloc(&d_res, sizeof(REAL));
//	hipMalloc(&d_a, n * sizeof(REAL));
//	hipMalloc(&d_b, n * sizeof(REAL));
//
//	hipMemset(d_res, 0, sizeof(REAL));
//	hipMemcpy(d_a, a, n * sizeof(REAL), hipMemcpyHostToDevice);
//	hipMemcpy(d_b, b, n * sizeof(REAL), hipMemcpyHostToDevice);
//
//	dotKernel << <(n + TPB - 1) / TPB, TPB >> >(d_res, d_a, d_b, n);
//	hipMemcpy(res, d_res, sizeof(REAL), hipMemcpyDeviceToHost);
//
//	hipFree(d_res);
//	hipFree(d_a);
//	hipFree(d_b);
//}


int divUp(int a, int b) { return (a + b - 1) / b; }

__device__
int idxClip(int idx, int idxMax) 
{
   return idx >(idxMax - 1) ? (idxMax - 1) : (idx < 0 ? 0 : idx);
   //return idx >(idxMax - 1) ? (idxMax - 1) : (idx );
}

__device__
int flatten(int col, int row, int width, int height) 
{
   return idxClip(col, width) + idxClip(row, height)*width;
}


/////// 第一种relax_p_kernel（）写法
__global__
void relax_p_kernel(REAL *d_pvec,  REAL *d_fvec, int w, int h, REAL ht)  // 测试通过！！！！
{
  HIP_DYNAMIC_SHARED( REAL, s_in)

  REAL ht2 = ht*ht;  // 多重网格中用于变化的空间步长

  // 设置扩充边界, device函数应该是简单的函数，计算中不能带循环，会特别影响速度！！！！
  // set_relax_p_bdry(d_bd_xl, d_bd_xr, d_bd_yl, d_bd_yr, d_pvec, w, h);

  // global indices
  int col = threadIdx.x + blockDim.x * blockIdx.x;
  int row = threadIdx.y + blockDim.y * blockIdx.y;

  //if ((col >= w) || (row >= h))
  //{
	 // printf("boundary error in relax_p_Kernel!!!\n ");
	 // return;
  //}

  //if ((col >= w) )
  //{
	 // printf("boundary error  col >= w \n ");
	 // printf("col = %d,  w = %d \n ", col, w);
	 // return;
  //}

  //if ((row >= h))
  //{
	 // printf("boundary error  row >= h \n ");
	 // printf("row = %d,  h = %d \n ", row, h);
	 // return;
  //}

  int idx = flatten(col, row, w, h);

  // local width and height
  int s_w = blockDim.x + 2 * RAD;
  int s_h = blockDim.y + 2 * RAD;

  // local indices
  int s_col = threadIdx.x + RAD;
  int s_row = threadIdx.y + RAD;
  int s_idx = flatten(s_col, s_row, s_w, s_h);
  
  // Load regular cells
  s_in[s_idx] = d_pvec[idx];

  // Load halo cells
  if (threadIdx.x < RAD) 
  {
    s_in[flatten(s_col - RAD, s_row, s_w, s_h)]        = d_pvec[flatten(col - RAD, row, w, h)];
    s_in[flatten(s_col + blockDim.x, s_row, s_w, s_h)] = d_pvec[flatten(col + blockDim.x, row, w, h)];
  }
  if (threadIdx.y < RAD) 
  {
    s_in[flatten(s_col, s_row - RAD, s_w, s_h)]        = d_pvec[flatten(col, row - RAD, w, h)];
    s_in[flatten(s_col, s_row + blockDim.y, s_w, s_h)] = d_pvec[flatten(col, row + blockDim.y, w, h)];
  }

  if (row == 0) 
  {
	  //int sid = flatten(col, 1, w, h);
	  //d_pvec[idx] = d_pvec[sid];

	  d_pvec[idx] = 0.0;

	  return;      // 注意: return 不能省略!!!
  }

  if (row == h - 1)     
  {  
	  //int sid = flatten(col, h-2, w, h);
	  //d_pvec[idx] = d_pvec[sid];

	  d_pvec[idx] = 0.0;

	  return;  // 注意: return 不能省略!!!
  }

  if (col == 0)
  {	  
	  //int sid = flatten(w-2, row, w, h);
	  //d_pvec[idx] = d_pvec[sid];

	  d_pvec[idx] = 0.0;

	  return;  // 注意: return 不能省略!!!
  }

  if ( col == w - 1)
  {
	  //int sid = flatten(1, row, w, h);
	  //d_pvec[idx] = d_pvec[sid];   

	  d_pvec[idx] = 0.0;

	  return;  // 注意: return 不能省略!!!
  }

  __syncthreads();

  // 迭代公式: u[i][j] = 0.25*(u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1]) - 0.25*(h*h)*f[i][j];
  REAL tp = 0.25 * (s_in[flatten(s_col - 1, s_row, s_w, s_h)]
                   + s_in[flatten(s_col + 1, s_row, s_w, s_h)] 
                   + s_in[flatten(s_col, s_row - 1, s_w, s_h)] 
                   + s_in[flatten(s_col, s_row + 1, s_w, s_h)]);
  d_pvec[idx] = tp - 0.25*d_fvec[idx] * ht2;

  //这种写法和上面前一种效果一样
  //REAL coef = -4.0 / ht2;
  //REAL src = -(s_in[flatten(s_col - 1, s_row, s_w, s_h)]
	 //         + s_in[flatten(s_col + 1, s_row, s_w, s_h)]
	 //         + s_in[flatten(s_col, s_row - 1, s_w, s_h)]
	 //         + s_in[flatten(s_col, s_row + 1, s_w, s_h)]) / ht2 
	 //         + d_fvec[flatten(col, row, w, h)]; 
  //d_pvec[idx] = src / coef;

  //这种写法和上面前两种效果一样
  //REAL coef = -4.0 / ht2;
  //REAL src = -(s_in[flatten(s_col - 1, s_row, s_w, s_h)]
	 //         + s_in[flatten(s_col + 1, s_row, s_w, s_h)]
	 //         + s_in[flatten(s_col, s_row - 1, s_w, s_h)]
	 //         + s_in[flatten(s_col, s_row + 1, s_w, s_h)]) / ht2;
  //d_pvec[idx] = (src+ d_fvec[idx]) / coef;

}



__global__
void relax_p_kernel_no_shared_memry(REAL *d_pvec, REAL *d_fvec, int w, int h, REAL ht) 
{
	//HIP_DYNAMIC_SHARED( REAL, s_in)

	REAL ht2 = ht*ht;  // 多重网格中用于变化的空间步长

	// 设置扩充边界, device函数应该是简单的函数，计算中不能带循环，会特别影响速度！！！！
	// set_relax_p_bdry(d_bd_xl, d_bd_xr, d_bd_yl, d_bd_yr, d_pvec, w, h);

	// global indices
	int col = threadIdx.x + blockDim.x * blockIdx.x;
	int row = threadIdx.y + blockDim.y * blockIdx.y;

	int idx = flatten(col, row, w, h);

	int idx0 = flatten(col - 1, row,     w, h);
	int idx1 = flatten(col,     row - 1, w, h);	
	int idx2 = flatten(col + 1, row,     w, h);
	int idx3 = flatten(col    , row + 1, w, h);

	if (row == 0)
	{
		//int sid = flatten(col, 1, w, h);
		//d_pvec[idx] = d_pvec[sid];

		d_pvec[idx] = 0.0;

		return;      // 注意: return 不能省略!!!
	}

	if (row == h - 1)
	{
		//int sid = flatten(col, h-2, w, h);
		//d_pvec[idx] = d_pvec[sid];

		d_pvec[idx] = 0.0;

		return;  // 注意: return 不能省略!!!
	}

	if (col == 0)
	{
		//int sid = flatten(w-2, row, w, h);
		//d_pvec[idx] = d_pvec[sid];

		d_pvec[idx] = 0.0;

		return;  // 注意: return 不能省略!!!
	}

	if (col == w - 1)
	{
		//int sid = flatten(1, row, w, h);
		//d_pvec[idx] = d_pvec[sid];   

		d_pvec[idx] = 0.0;

		return;  // 注意: return 不能省略!!!
	}

	__syncthreads();

	// //迭代公式: u[i][j] = 0.25*(u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1]) - 0.25*(h*h)*f[i][j];
	//REAL tp = 0.25 * (d_pvec[idx0] + d_pvec[idx1] + d_pvec[idx2] + d_pvec[idx3]);

	//d_pvec[idx] = tp - 0.25*d_fvec[idx] * ht2;


	////////////////////liu//////////////////
	if (row < h - 1 && col < w - 1 && row>0 && col>0 && (row + col) % 2 == 0)
	{
		d_pvec[idx] = 0.25 * (d_pvec[idx0] + d_pvec[idx1] + d_pvec[idx2] + d_pvec[idx3]- d_fvec[idx] * ht2);
	}
	__syncthreads();
	if (row < h - 1 && col < w - 1 && row>0 && col>0 && (row + col) % 2 == 1)
	{
		d_pvec[idx] = 0.25 * (d_pvec[idx0] + d_pvec[idx1] + d_pvec[idx2] + d_pvec[idx3] - d_fvec[idx] * ht2);
	}
	__syncthreads();

	///////////////////////////////////////////////////////////

	//if (idx < 256)
	//{
	//	printf("idx=%d  d_fvec[%d]=%lf\n  ",idx,idx,d_fvec[idx]);
	//}
	

	//这种写法和上面前一种效果一样
	//REAL coef = -4.0 / ht2;
	//REAL src = -(s_in[flatten(s_col - 1, s_row, s_w, s_h)]
	//         + s_in[flatten(s_col + 1, s_row, s_w, s_h)]
	//         + s_in[flatten(s_col, s_row - 1, s_w, s_h)]
	//         + s_in[flatten(s_col, s_row + 1, s_w, s_h)]) / ht2 
	//         + d_fvec[flatten(col, row, w, h)]; 
	//d_pvec[idx] = src / coef;

	//这种写法和上面前两种效果一样
	//REAL coef = -4.0 / ht2;
	//REAL src = -(s_in[flatten(s_col - 1, s_row, s_w, s_h)]
	//         + s_in[flatten(s_col + 1, s_row, s_w, s_h)]
	//         + s_in[flatten(s_col, s_row - 1, s_w, s_h)]
	//         + s_in[flatten(s_col, s_row + 1, s_w, s_h)]) / ht2;
	//d_pvec[idx] = (src+ d_fvec[idx]) / coef;

}


// 完成 2021.07.08, 有问题, debugging........  
//__global__
//void relax_p_kernel(REAL *d_p, REAL *d_f, int w, int h, int gnx)  // 输入输出数组为d_p, 计算完成后会更新值，d_p = [0...nxt-1, 0...nyt-1], 
//														            // 输入数组为d_f， d_f = [0...nxt-1, 0...nyt-1],
//														            // w和h给定的是数组实际维数，d_p和d_f是nxt*nyt维的矩阵。(col,row)二维索引从(0,0)开始。     
//{
//
//	int i, j, iter;
//	REAL  coef, src;
//
//    REAL ht = (REAL)gnx / w;  // 多重网格中用于变化的空间步长
//	REAL ht2 = ht*ht;
//
//	HIP_DYNAMIC_SHARED( REAL, s_in)
//
//	// global indices
//	const int col = threadIdx.x + blockDim.x * blockIdx.x;
//	const int row = threadIdx.y + blockDim.y * blockIdx.y;
//
//	if ((col >= w) || (row >= h)) { printf("boundary errow for col and row\n");  return;}
//
//	 int idx = flatten(col, row, w, h);
//	// local width and height
//
//	const int s_w = blockDim.x + 2 * RAD;
//	const int s_h = blockDim.y + 2 * RAD;
//
//	// local indices
//	const int s_col = threadIdx.x + RAD;
//	const int s_row = threadIdx.y + RAD;
//	const int s_idx = flatten(s_col, s_row, s_w, s_h);
//
//	// Load regular cells
//	s_in[s_idx] = d_p[idx];
//
//	// Load halo cells
//	if (threadIdx.x < RAD)
//	{
//		int local_flat_xm = flatten(s_col - RAD,        s_row, s_w, s_h);
//		int local_flat_xp = flatten(s_col + blockDim.x, s_row, s_w, s_h);
//		int       flat_xm = flatten(col - RAD, row,         w,        h);
//		int       flat_xp = flatten(col + blockDim.x,     row,   w,   h);
//
//		s_in[local_flat_xm] = d_p[flat_xm];
//		s_in[local_flat_xp] = d_p[flat_xp];
//	}
//
//	if (threadIdx.y < RAD)
//	{
//		int local_flat_ym = flatten(s_col, s_row - RAD,        s_w, s_h);
//		int local_flat_yp = flatten(s_col, s_row + blockDim.y, s_w, s_h);
//		int       flat_ym = flatten(col, row - RAD,              w,   h);
//		int       flat_yp = flatten(col, row + blockDim.y,       w,   h);
//
//		s_in[local_flat_ym] = d_p[flat_ym];
//		s_in[local_flat_yp] = d_p[flat_yp];
//	}
//
//	src = d_f[idx];    // 需要检查确定变量src能不能累加!!!
//
//	printf("d_f[%d] = %f\n", idx, d_f[idx]);
//
//	if (col > 0 && col < w - 1)
//	{
//		int xm1 = flatten(s_col - 1, s_row, s_w, s_h);
//		int xp1 = flatten(s_col + 1, s_row, s_w, s_h);
//
//		src -= (s_in[xm1] + s_in[xp1]) / ht2;
//		coef = 2.0 / ht2;
//
//		return;
//	}
//	else if (col == 0)
//	{
//		int id0 = flatten(1, row, w, h);
//		int id1 = flatten(w - 1, row, w, h);
//
//		src -= (d_p[id0] + d_p[id1]) / ht2;
//		coef = 2.0 / ht2;
//
//		return;
//	}
//	else if (col == w - 1)
//	{
//		int id0 = flatten(0, row, w, h);
//		int id1 = flatten(w - 2, row, w, h);
//
//		src -= (d_p[id0] + d_p[id1]) / ht2;
//		coef = 2.0 / ht2;
//
//		return;
//	}
//
//
//	if (row > 0 && row < h - 1)
//	{
//		int ym1 = flatten(s_col, s_row - 1, s_w, s_h);
//		int yp1 = flatten(s_col, s_row + 1, s_w, s_h);
//
//		src -= (s_in[ym1] + s_in[yp1]) / ht2;
//		coef += 2.0 / ht2;
//
//		return;
//	}
//	else if (row == 0)
//	{
//		int id = flatten(col, 0, w, h);
//		src -= d_p[id] / ht2;
//		coef += 1.0 / ht2;
//	}
//	else if (row == h - 1)
//	{
//		int id = flatten(col, h - 2, w, h);
//		src -= d_p[id] / ht2;
//		coef += 1.0 / ht2;
//
//		return;
//	}
//
//	__syncthreads();
//
//	d_p[idx] = -src / coef;
//
//}



/**************************** begin  已修改的 *************************************/


void kernelLauncher2(REAL *d_pvec, REAL  *d_fvec, int w, int h)  // 还有问题, 原子操作atomicMax(),cuda不识别!  测试 mat_max_kernel()的包装函数!
{
	const dim3 blockSize(TX, TY);
	const dim3 gridSize(divUp(w, TX), divUp(h, TY));
	//const size_t smSz = (TX + 2 * RAD)*(TY + 2 * RAD) * sizeof(REAL);

	//mat_max_kernel << <gridSize, blockSize >> >(d_pvec, d_fvec, w, h);
}

//////   test kernel !!!!!
//void test_kernelLauncher(REAL *d_pvec, REAL  *d_fvec, int w, int h, REAL ht)
//{
//	const dim3 blockSize(TX, TY);
//	const dim3 gridSize(divUp(w, TX), divUp(h, TY));
//	const size_t smSz = (TX + 2 * RAD)*(TY + 2 * RAD) * sizeof(REAL);
//
//	//restrict_kernel << <gridSize, blockSize >> >(d_pvec, d_fvec, w, h);  // 测试通过
//
//	//prolong_kernel << <gridSize, blockSize >> >(d_fvec, d_pvec, w, h);   // 测试通过
//
//	// zero_matrix_kernel << <gridSize, blockSize>> >(d_pvec, w, h);   // 测试通过
//
//	// hipMemset(d_pvec, 0, (w *h ) * sizeof(REAL)); // 测试通过
//
//	//mat_add_kernel << <gridSize, blockSize >> > (d_pvec, d_fvec, d_pvec, w, h);  // 测试通过 ！
//
//	// mat_sub_kernel << <gridSize, blockSize >> > (d_pvec, d_fvec, d_pvec, w, h); // 测试通过 ！
//
//	// hipLaunchKernelGGL(laplace_p_kernel, dim3(gridSize), dim3(blockSize), smSz , 0, d_pvec, d_fvec, w, h, nx);  // 测试通过, 注意这个求拉普拉斯p的边界条件!
//
//	// mat_copy_kernel << <gridSize, blockSize >> > (d_pvec, d_fvec, w, h); // 测试通过 !
//
//	//relax_p_kernel << <gridSize, blockSize, smSz >> >(d_pvec, d_fvec, w, h, ht);
//}


void relax_p_kernelLauncher(REAL *d_pvec, REAL  *d_fvec, int w, int h, REAL ht)
{
	const dim3   blockSize(TX, TY);
	const dim3   gridSize(divUp(w, TX), divUp(h, TY));
	const size_t smSz = (TX + 2 * RAD)*(TY + 2 * RAD) * sizeof(REAL);

	//// 开辟对应扩充边界的一维显存
	//REAL *d_bd_xl, *d_bd_xr, *d_bd_yl, *d_bd_yr;
	//hipMalloc(&d_bd_xl, w * sizeof(REAL));
	//hipMalloc(&d_bd_xr, w * sizeof(REAL));
	//hipMalloc(&d_bd_yl, h * sizeof(REAL));
	//hipMalloc(&d_bd_yr, h * sizeof(REAL));

	// 核心计算
	//relax_p_kernel << <gridSize, blockSize, smSz >> >(d_pvec, d_fvec, w, h, ht);
	hipLaunchKernelGGL(relax_p_kernel_no_shared_memry,gridSize,blockSize,0,0,d_pvec,d_fvec,w,h,ht);
	//relax_p_kernel_no_shared_memry << <gridSize, blockSize>> > (d_pvec, d_fvec, w, h, ht);


	//// 释放显存
	//hipFree(d_bd_xl);
	//hipFree(d_bd_xr);
	//hipFree(d_bd_yl);
	//hipFree(d_bd_yr);
}

__device__
void set_relax_p_bdry(REAL *bd_xl, REAL *bd_xr, REAL *bd_yl, REAL *bd_yr, REAL *d_pvec, int w, int h)
{
	// 设置扩展边界
	for (int i = 0; i < w; i++)
	{
		//bd_xl[i] = pmatf[i][1];
		//bd_xr[i] = pmatf[i][nyf];
		bd_xl[i] = d_pvec[flatten(i, 0 ,  w, h)];
		bd_xr[i] = d_pvec[flatten(i, h-1, w, h)];
	}

	for (int j = 0; j < h; j++)
	{
		//bd_yu[j] = pmatf[nxf][j];
		//bd_yd[j] = pmatf[1][j];
		bd_yl[j] = d_pvec[flatten(w-1 , j, w, h)];
		bd_yr[j] = d_pvec[flatten(0,    j, w, h)];
	}

}


void test_restrict_Launcher(REAL *d_ufine, REAL  *d_ucorse, int w, int h)
{
	const dim3 blockSize(TX, TY);
	const dim3 gridSize(divUp(w, TX), divUp(h, TY));
	//const size_t smSz = (TX + 2 * RAD)*(TY + 2 * RAD) * sizeof(REAL);

	hipLaunchKernelGGL(restrict_kernel,gridSize,blockSize,0,0,d_ufine,d_ucorse,w,h);
	//restrict_kernel << <gridSize, blockSize >> >(d_ufine, d_ucorse, w, h); // 测试通过
}


__global__
void restrict_kernel(REAL *d_ufine, REAL *d_ucorse, int w, int h)  // d_ucorse[0...w*h-1], d_ufine[0...(2w)*(2h)-1]   // 测试通过
{

	// global indices
	const int col = threadIdx.x + blockDim.x * blockIdx.x;
	const int row = threadIdx.y + blockDim.y * blockIdx.y;

	if ((col >= w) || (row >= h))
	{
		printf("boundary error in restrict_kernel()!!!\n ");
		return;
	}

	const int idx = flatten(col, row, w, h);

	//printf("p[%d] = %f\n", idx, d_ufine[idx]);

	const int id0 = flatten(2 * col, 2 * row, 2 * w, 2 * h);
	const int id1 = flatten(2 * col + 1, 2 * row, 2 * w, 2 * h);
	const int id2 = flatten(2 * col, 2 * row + 1, 2 * w, 2 * h);
	const int id3 = flatten(2 * col + 1, 2 * row + 1, 2 * w, 2 * h);

	//syncthreads();

	d_ucorse[idx] = 0.25 * (d_ufine[id0] + d_ufine[id1] + d_ufine[id2] + d_ufine[id3]);
}


__global__
void prolong_kernel(REAL *d_ucorse, REAL *d_ufine, int w, int h) // d_ucorse[0...w*h-1], d_ufine[0...(2w)*(2h)-1]   // 测试通过
{
	// global indices
	const int col = threadIdx.x + blockDim.x * blockIdx.x;
	const int row = threadIdx.y + blockDim.y * blockIdx.y;

	const int idx = flatten(col, row, w, h);

	const int id0 = flatten(2 * col + 1, 2 * row + 1, 2 * w, 2 * h);
	const int id1 = flatten(2 * col + 1, 2 * row, 2 * w, 2 * h);
	const int id2 = flatten(2 * col, 2 * row + 1, 2 * w, 2 * h);
	const int id3 = flatten(2 * col, 2 * row, 2 * w, 2 * h);

	d_ufine[id0] = d_ufine[id1] = d_ufine[id2] = d_ufine[id3] = d_ucorse[idx];

	return;
}


//// 用法: REAL *d_max;    // 还有问题, 原子操作atomicMax(),cuda不识别!   
//__global__
//void mat_max_kernel(REAL *d_max, REAL *d_a, int w, int h) // 找出数组d_a中绝对值最大的数,并返回这个最大的绝对值在指针d_max中   // 需要做测试 ！！！！！！！！！！！！！！！！！！！！！！！！！！！！
//{
//	// global indices
//	const int col = threadIdx.x + blockDim.x * blockIdx.x;
//	const int row = threadIdx.y + blockDim.y * blockIdx.y;
//
//	if ((col >= w) || (row >= h))
//	{
//		printf("boundary error in mat_max_kernel()!!!\n ");
//		return;
//	}
//
//	extern __shared__ REAL ss_in[TX*TY];
//
//	const int idx = flatten(col, row, w, h);
//	
//	const int s_w = blockDim.x;
//	const int s_h = blockDim.y;
//
//	const int s_col = threadIdx.x;
//	const int s_row = threadIdx.y;
//
//	const int s_idx = flatten(s_col, s_row, s_w, s_h);
//
//	ss_in[s_idx] = d_a[idx];
//
//	REAL x = 0.0;
//
//	for (int j = 0; j < blockDim.x * blockDim.y; j++)
//	{
//		if (ss_in[j] > x) x = ss_in[j];
//		else if ((-ss_in[j]) > x)  x = -ss_in[j];
//		//return x;
//	}
//
//	atomicMax(d_max, x);
//}


__global__
void mat_copy_kernel(REAL *d_a, REAL *d_b, int w, int h)    // 测试通过！！！！注意这个拷贝函数用于两个显存数组之间的拷贝!
                                                                // 不同于hipMemcpy()函数,hipMemcpy()用于主存与显存之间的拷贝! 
{
	// global indices
	const int col = threadIdx.x + blockDim.x * blockIdx.x;
	const int row = threadIdx.y + blockDim.y * blockIdx.y;
	
	if ((col >= w) || (row >= h))
	{
		printf("boundary error in mat_copy_kernel()!!!\n ");
		return;
	}

	const int idx = flatten(col, row, w, h);

	d_a[idx] = d_b[idx];

	return;
}


__global__
void laplace_p_kernel(REAL *d_pvec, REAL  *d_lap_p, int w, int h, int nx)   // // 测试通过, 注意这个求拉普拉斯p的边界条件!
{
	HIP_DYNAMIC_SHARED( REAL, slp_in)

	REAL ht = (REAL) nx / w;   // 多重网格中用于变化的空间步长
	REAL ht2 = ht*ht;

	// global indices
	const int col = threadIdx.x + blockDim.x * blockIdx.x;
	const int row = threadIdx.y + blockDim.y * blockIdx.y;

	if ((col >= w) || (row >= h))
	{
		printf("boundary error in laplace_p_kernel()\n ");
		return;
	}

	const int idx = flatten(col, row, w, h);

	// local width and height
	const int s_w = blockDim.x + 2 * RAD;
	const int s_h = blockDim.y + 2 * RAD;

	// local indices
	const int s_col = threadIdx.x + RAD;
	const int s_row = threadIdx.y + RAD;
	const int s_idx = flatten(s_col, s_row, s_w, s_h);

	// Load regular cells
	slp_in[s_idx] = d_pvec[idx];

	// Load halo cells
	if (threadIdx.x < RAD)
	{
		slp_in[flatten(s_col - RAD, s_row, s_w, s_h)]        = d_pvec[flatten(col - RAD, row, w, h)];
		slp_in[flatten(s_col + blockDim.x, s_row, s_w, s_h)] = d_pvec[flatten(col + blockDim.x, row, w, h)];
	}
	if (threadIdx.y < RAD)
	{
		slp_in[flatten(s_col, s_row - RAD, s_w, s_h)]        = d_pvec[flatten(col, row - RAD, w, h)];
		slp_in[flatten(s_col, s_row + blockDim.y, s_w, s_h)] = d_pvec[flatten(col, row + blockDim.y, w, h)];
	}

	if ((col == 0))
	{
		d_lap_p[idx] = 0;  // 赋固定值, 因为是纽曼边界条件
		return;      // 注意: return 不能省略!!!
	}

	if (col == w - 1)
	{
		d_lap_p[idx] = 0;   // 赋固定值, 因为是纽曼边界条件
		return;  // 注意: return 不能省略!!!
	}

	if (row == 0)      // 周期边界
	{
		int sid0 = flatten(col, 1,     w, h);
		int sid1 = flatten(col, h - 2, w, h);
		int sid2 = flatten(col, 0,     w, h);
		
		d_lap_p[idx] = (d_pvec[sid0] + d_pvec[sid1] - 2.0*d_pvec[sid2])/ht2;  

		return;  // 注意: return 不能省略!!!
	}

	if (row == h - 1)   // 周期边界
	{
		int sid0 = flatten(col, 0, w, h);
		int sid1 = flatten(col, h - 2, w, h);
		int sid2 = flatten(col, h-1, w, h);

		d_lap_p[idx] = (d_pvec[sid0] + d_pvec[sid1] - 2.0 * d_pvec[sid2]) / ht2;
		return;                          // 注意: return 不能省略!!!
	}

	__syncthreads();

	// Laplace公式 lap[i][j] = (u[i+][j] + u[i-1][j] + u[i][j+1] + u[i][j-1] - 4.0* u[i][j])/(h*h)
	REAL tp = (slp_in[flatten(s_col - 1, s_row, s_w, s_h)]
		     +  slp_in[flatten(s_col + 1, s_row, s_w, s_h)]
		     +  slp_in[flatten(s_col, s_row - 1, s_w, s_h)]
		     +  slp_in[flatten(s_col, s_row + 1, s_w, s_h)]
		     -  4.0*slp_in[flatten(s_col, s_row, s_w, s_h)]) / ht2;

	d_lap_p[idx] = tp;

}


// 调用方法: mat_sub_kernel << <gridSize, blockSize, smSz >> >(d_a, d_b, d_c, w, h);
__global__
void mat_sub_kernel(REAL *d_a, REAL *d_b, REAL *d_c, int w, int h)   // 测试通过 ！
{
	// global indices
	const int col = threadIdx.x + blockDim.x * blockIdx.x;
	const int row = threadIdx.y + blockDim.y * blockIdx.y;
	
	if ((col >= w) || (row >= h))
	{
		printf("boundary error in mat_sub_kernel()!!!\n ");
		return;
	}

	const int idx = flatten(col, row, w, h);

	d_a[idx] = d_b[idx] - d_c[idx];

	return;
}


// 调用方法: mat_add_kernel << <gridSize, blockSize, smSz >> >(d_a, d_b, d_c, w, h);
__global__
void mat_add_kernel(REAL *d_a, REAL *d_b, REAL *d_c, int w, int h)      // 测试通过
{
	// global indices
	const int col = threadIdx.x + blockDim.x * blockIdx.x;
	const int row = threadIdx.y + blockDim.y * blockIdx.y;
	
	if ((col >= w) || (row >= h))
	{
		printf("boundary error in mat_add_kernel()!!!\n ");
		return;
	}

	const int idx = flatten(col, row, w, h);

	d_a[idx] = d_b[idx] + d_c[idx];

	return;
}


// 应该也可以直接使用  hipMemset(d_a, 0, N*sizeof(REAL)); 
// 调用方法:zero_matrix_kernel << <gridSize, blockSize, smSz >> >(d_a, w, h);    // 测试通过
__global__
void zero_matrix_kernel(REAL *d_a, int w, int h)
{
	// global indices
	const int col = threadIdx.x + blockDim.x * blockIdx.x;
	const int row = threadIdx.y + blockDim.y * blockIdx.y;

	if ((col >= w) || (row >= h))
	{
		printf("boundary error in zero_matrix_kernel()!!!\n ");
		return;
	}

	const int idx = flatten(col, row, w, h);

	d_a[idx] = 0.0;

	return;
}



/******************* end 已修改的 ************************/



/**********  还没修改
void vcycle_uv(REAL **uf, REAL **ff, int nxf, int nyf, int ilevel)
{
	extern int n_level;

	relax_p(uf, ff, nxf, nyf);


	if (ilevel < n_level)
	{
		int nxc, nyc;
		REAL **rf, **uc, **fc;

		nxc = nxf / 2;
		nyc = nyf / 2;

		rf = dmatrix(1, nxf, 1, nyf);
		uc = dmatrix(1, nxc, 1, nyc);
		fc = dmatrix(1, nxc, 1, nyc);

		residual_p(rf, uf, ff, nxf, nyf);

		restrict1(rf, fc, nxc, nyc);
		zero_matrix(uc, 1, nxc, 1, nyc);
		vcycle_uv(uc, fc, nxc, nyc, ilevel + 1);

		prolong(uc, rf, nxc, nyc);
		mat_add(uf, uf, rf, 1, nxf, 1, nyf);

		relax_p(uf, ff, nxf, nyf);

		free_dmatrix(rf, 1, nxf, 1, nyf);
		free_dmatrix(uc, 1, nxc, 1, nyc);
		free_dmatrix(fc, 1, nxc, 1, nyc);
	}
}

void MG_Poisson(REAL **p, REAL **f)     // 多重网格方法求解压力泊松方程
{
	extern int nx, ny;
	extern REAL **sor, **workv;

	int    max_it = 2000, it_mg = 1;
	REAL tol = 1.0e-5, resid = 1.0;

	mat_copy(workv, p, 1, nx, 1, ny);

	while (it_mg <= max_it && resid >= tol)
	{

		vcycle_uv(p, f, nx, ny, 1);

		pressure_update(p);

		ijloop
		{
			sor[i][j] = workv[i][j] - p[i][j];  // mat_sub
		}

		resid = mat_max(sor, 1, nx, 1, ny);
		mat_copy(workv, p, 1, nx, 1, ny);

		it_mg++;
	}

	printf("Mac iteration = %d  residual = %16.15f \n", it_mg, resid);

	return;
}


void residual_p(REAL **r, REAL **u, REAL **f, int nxt, int nyt)
{
	laplace_p(u, r, nxt, nyt);
	mat_sub(r, f, r, 1, nxt, 1, nyt);
}


void laplace_p(REAL **p, REAL **lap_p, int nxt, int nyt)
{
	REAL **dpdx, **dpdy;

	dpdx = dmatrix(0, nxt, 1, nyt);
	dpdy = dmatrix(1, nxt, 0, nyt);

	grad_p(p, dpdx, dpdy, nxt, nyt);
	div_uv(dpdx, dpdy, lap_p, nxt, nyt);

	free_dmatrix(dpdx, 0, nxt, 1, nyt);
	free_dmatrix(dpdy, 1, nxt, 0, nyt);

}

********/
